用TTree 结构存储事件数据
一、学习目的¶
学习如何将事件以TTree的格式存储到ROOT文件
二、学习内容¶
束流入射到靶上产生核反应，靶上生成的中子和$\gamma$入射到长条型塑料闪烁体探测器的不同位置，并在不同深度与探测介质产生相互作用。
束流的时间信息$T_{Beam}$由Timing Scintillator(薄塑料闪烁体+PMT)给出。
入射粒子形成的光信号向探测器两端传输，经PMT放大形成输出信号。
输出信号分成两路，一路经CFD甄别后形成时间信号 tu 和 td ,送入TDC的stop端；另一路延迟送入QDC，得到能量信号 Qu,Qd 。
获取系统的触发为 $T_u\times T_d\times T_{Beam}$ , 触发信号中$T_{Beam}$的时间为主时序。
触发信号经适当调节送到TDC的Start端和QDC的Gate端。
将入射粒子的原始信息，以及探测器信息逐事件存入root文件的TTree结构中.
入射粒子原始信息：粒子种类、能量，粒子在探测器上的入射位置、入射深度等
探测器信息: 两端时间信息 tu 、 td ，能量沉积信息 Qu , Qd
打开root文件，进行基本的分析。
三、需要处理的问题——中子飞行时间谱测量原理¶
束流¶
探测器参数¶
闪烁体探测器长度：2L
光在闪烁体的传播速度: $v_{sc}$，光衰减常数：$\lambda$
$d=\sqrt{D^2+x^2}$
中子飞行时间：$TOF_n(ns)=\frac{72 \cdot d(m)}{\sqrt{E_n(MeV)}}$, gamma飞行时间：$TOF_{\gamma}(ns)=3\cdot d(m)$
探测器两端时间：$t_u=TOF+(L-x)/v_{sc}$， $t_d=TOF+(L+x)/v_{sc}$
探测器两端能量(x 处能量沉积Q)：$Q_u=Q_0 exp[-(L-x)/\lambda]$， $Q_d=Q_0 exp[-(L+x)/\lambda]$
物理过程推导¶
四、具体的代码实施过程¶
具体步骤¶
1. 声明所需的变量；
2. 建立TTree的存储；
3. 建立具体的分支；
4. 单个时间测存储；
5. 循环时间的存储；
6. 修改为root脚本运行。声明所需变量¶
In [1]:
//声明需要注意的问题：
//1、“同样可以使用const double aa=500.；进行声明”，其中加分号终端不会输出；
//2、100后不加点为整型输出，加小数点为浮点型。
//3、可以在后面一直加声明，不需要再阐述数据类型，但是不建议太多，容易乱
const Double_t Dis=500.,Len=100.,del_Dis=5.; //单位为cm，
const Double_t Lambda=380.; //cm, 光的衰减常数
const Double_t V_sc=7.5; //ns/cm, 光在闪烁体中的传播速度
const Double_t Rg=0.3; //出射gamma射线和中子的比值 //假设为均匀分布
const Double_t En0=100.; //MeV; 中子初始能量
const Double_t Eg0=1; //MeV; γ射线的能量
const Double_t EnRes=50.;//MeV
const Double_t TiRes=1.; //ns
const Double_t QRes=0.1;//relative energy resolution(FWHM) of the scin.
Double_t Tu_off=5.5; //time offset，//束流飞行时间+PMT的渡越时间+电缆上的传输时间
Double_t Td_off=20.4; //time offset
In [2]:
//声明所需要的参数，即需要输出的参数
//* 入射粒子原始信息：粒子种类、能量，粒子在探测器上的入射位置、入射深度等
//* 探测器信息: 两端时间信息 tu 、 td ，能量沉积信息 Qu , Qd
Double_t X,Energy; //应该是均匀分布
Int_t PID; //pid=1 is neutron、pid=0 is gamma
Double_t T_u,T_d,Q_u,Q_d;
Double_t TOF,cTOF; //TOF:粒子实际飞行时间，cTOF：计算得到的TOF
In [3]:
//声明中间参量，不需要输出
Double_t Deepth; //粒子的飞行深度
Double_t Dis_Deepth;
Double_t Range;
定义TTree¶
In [4]:
//新建tree的文件，然后声明tree
TFile *Output_File=new TFile("tree.root","recreate");
TTree *Output_Tree=new TTree("tree","tree structure");
In [5]:
Output_Tree->Branch("PID",&PID,"PID/I");
Output_Tree->Branch("X",&X,"X/D");
Output_Tree->Branch("Energy",&Energy,"Energy/D");
Output_Tree->Branch("T_u",&T_u,"T_u/D");
Output_Tree->Branch("T_d",&T_d,"T_d/D");
Output_Tree->Branch("Q_u",&Q_u,"Q_u/D");
Output_Tree->Branch("Q_d",&Q_d,"Q_d/D");
Output_Tree->Branch("TOF",&TOF,"TOF/D");
Output_Tree->Branch("cTOF",&cTOF,"cTOF/D");
新建一个直方图，用来存储飞行时间¶
In [6]:
TH1D *hctof=new TH1D("hctof","neutron time of flight",1000,0,100);
声明随机数¶
In [7]:
TRandom3 *gr=new TRandom3(0);
In [8]:
//开启循环
for(Int_t i_event=0;i_event<100000;i_event++){
X=gr->Uniform(-Len,Len);
Deepth=gr->Uniform(-del_Dis/2,del_Dis/2);
Dis_Deepth=Dis+Deepth;
Range=TMath::Sqrt(X*X+Dis_Deepth*Dis_Deepth);
//判断为哪种粒子
if(gr->Uniform()<Rg){
PID=0;
Energy=Eg0;
TOF=10/3*(Range*0.01);
}
else{
PID=1;
Energy=gr->Gaus(En0,EnRes/2.355);
TOF=72.29824/TMath::Sqrt(Energy)*(Range*0.01);
}
//该给出时间信号了
T_u=TOF+(Len-X)/V_sc+gr->Gaus(0,TiRes/2.35)+Tu_off;
T_d=TOF+(Len+X)/V_sc+gr->Gaus(0,TiRes/2.35)+Td_off;
cTOF=(T_u+T_d)/2.;//simplified calculation.
hctof->Fill(cTOF);
//该存储能量信号了
Double_t q0=Energy*gr->Uniform();
Q_u=q0*TMath::Exp(-(Len-X)/Lambda);
Q_u=gr->Gaus(Q_u,Q_u*QRes/2.35);//QRes, relative energy resolution
Q_d=q0*TMath::Exp(-(Len+X)/Lambda);
Q_d=gr->Gaus(Q_d,Q_d*QRes/2.35);
Output_Tree->Fill();
}
hctof->Write();//写入预定义的histogram到文件
Output_Tree->Write();//写入TTree到文件
Output_File->Close();//关闭文件
观察是否一致¶
In [9]:
%jsroot on
In [10]:
TFile *Input_File=new TFile("tree.root");
Input_File->ls()
In [11]:
tree->Print()
In [12]:
tree->Show(0);
In [13]:
tree->Show(3)
In [14]:
TCanvas *c1=new TCanvas();//* 在ROOT环境下可省略
c1->Clear();//* 在ROOT环境下可省略
tree->Draw("T_d-T_u>>htx(500,-20,50)");//位置一维分布
c1->Draw();//* 在ROOT环境下可省略
In [15]:
tree->Draw("TOF>>hh(1000,0,100)");//实际飞行时间
c1->SetLogy();
c1->Draw();
In [16]:
TH1D *hh=(TH1D*) Input_File->Get("hctof"); //在代码中得到文件内数据指针的方法，在ROOT环境下直接用 hctof->Draw()
hh->Draw();//显示文件内histogram
c1->Draw();
In [17]:
c1->SetLogy(0);
gStyle->SetPalette(1);
tree->Draw("cTOF:T_d-T_u>>(2000,-20,50,1000,0,100)","","colz");//二维图显示
c1->Draw()
In [19]:
tree->Draw("cTOF:T_d-T_u>>(2000,-20,50,50,40,45)","PID==0","colz");//加入条件,gamma
c1->Draw();//为什么在两侧上弯？
In [20]:
tree->Draw("T_d-T_u:X","","colz");//观察两个参量的关联
c1->Draw();
In [23]:
tree->Draw("log(Q_u/Q_d):X","","colz");
c1->Draw()
In [18]:
!jupyter nbconvert 1.1_Creat_TTree.ipynb --to html